This knitr document contains the code used to analyze and visualize data from the VIS camera system and water data. Knit-ing this document will generate an HTML document that contains both the embedded R code and the output.
Attach the required R packages
#library(qvalue, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)
## Warning: package 'ggplot2' was built under R version 3.1.3
library(lubridate, warn.conflicts = FALSE)
library(MASS, warn.conflicts = FALSE)
library(car, warn.conflicts = FALSE)
## Warning: package 'car' was built under R version 3.1.2
library(nlme, warn.conflicts = FALSE)
library(mvtnorm, warn.conflicts = FALSE)
library(grid, warn.conflicts = FALSE)
Download the complete VIS data set. There were a total of 6,399 snapshots with VIS image data, but the download only includes the 6,207 snapshots that were successfully processed by PlantCV. Failed snapshots generally are those that lack a plant (empty pot controls, dead plants, etc.)
#if (!file.exists('vis_snapshots.txt')) {
# download.file('http://files.figshare.com/1845362/vis_snapshots.txt',
# 'vis_snapshots.txt')
#}
Read data and format for analysis
vis.data = read.table(file="vis_snapshots_nocorrect.csv", sep=",", header=TRUE)
Planting date
planting_date = as.POSIXct("2013-11-26")
Add water treatment column coded in barcodes.
vis.data$treatment <- NA
vis.data$treatment[grep("AA", vis.data$plant_id)] <- 100
vis.data$treatment[grep("AB", vis.data$plant_id)] <- 0
vis.data$treatment[grep("AC", vis.data$plant_id)] <- 16
vis.data$treatment[grep("AD", vis.data$plant_id)] <- 33
vis.data$treatment[grep("AE", vis.data$plant_id)] <- 66
Add plant genotype column coded in barcodes.
vis.data$genotype <- NA
vis.data$genotype[grep("p1", vis.data$plant_id)] <- 'A10'
vis.data$genotype[grep("p2", vis.data$plant_id)] <- 'B100'
vis.data$genotype[grep("r1", vis.data$plant_id)] <- 'R20'
vis.data$genotype[grep("r2", vis.data$plant_id)] <- 'R70'
vis.data$genotype[grep("r3", vis.data$plant_id)] <- 'R98'
vis.data$genotype[grep("r4", vis.data$plant_id)] <- 'R102'
vis.data$genotype[grep("r5", vis.data$plant_id)] <- 'R128'
vis.data$genotype[grep("r6", vis.data$plant_id)] <- 'R133'
vis.data$genotype[grep("r7", vis.data$plant_id)] <- 'R161'
vis.data$genotype[grep("r8", vis.data$plant_id)] <- 'R187'
Add genotype x treatment group column
vis.data$group = paste(vis.data$genotype,'-',vis.data$treatment,sep='')
Add calendar-time data column using the Unix-time data
vis.data$date = as.POSIXct(vis.data$datetime, origin = "1970-01-01")
Calculate days after planting from planting data
vis.data$dap = as.numeric(vis.data$date - planting_date)
Use linear regression to create a simple model for using water volume to predict soil volume water content.
Download data for soil wet/dry weight and volume water content measurements.
if (!file.exists('soil_weigth_vwc.txt')) {
download.file('http://files.figshare.com/1939954/soil_weigth_vwc.txt',
'soil_weigth_vwc.txt')
}
Read the soil volume water content data
vwc.data = read.table(file="soil_weigth_vwc.txt", sep="\t", header=TRUE)
Calculate the average soil dry weight
mean(vwc.data$weight_dry)
## [1] 72.92138
Create a linear model for water volume and volume water content. Water volumes >= 260 appear to have saturated the soil water carrying capacity, so remove them from the model.
vwc.lm = lm(vwc_wet ~ water_vol, data=vwc.data[vwc.data$water_vol < 260,])
summary(vwc.lm)
##
## Call:
## lm(formula = vwc_wet ~ water_vol, data = vwc.data[vwc.data$water_vol <
## 260, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5426 -1.6196 0.3976 1.4214 4.0771
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.032863 1.389638 -2.182 0.0391 *
## water_vol 0.233197 0.008145 28.630 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.082 on 24 degrees of freedom
## Multiple R-squared: 0.9716, Adjusted R-squared: 0.9704
## F-statistic: 819.7 on 1 and 24 DF, p-value: < 2.2e-16
Plot the model
vwc.plot = ggplot(vwc.data[vwc.data$water_vol < 260,], aes(x=water_vol, y=vwc_wet)) +
geom_point(size=2) +
geom_smooth(method='lm', formula=y~x) +
scale_x_continuous("Water volume (ml)") +
scale_y_continuous("Volume water content (%)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
print(vwc.plot)
Predict volume water contents for the water treatment groups
treatment.water = data.frame(water_vol=c(217, 144.5, 72))
treatment.water$vwc = predict(object = vwc.lm, newdata=treatment.water)
print(treatment.water)
## water_vol vwc
## 1 217.0 47.57084
## 2 144.5 30.66407
## 3 72.0 13.75731
LemnaTec VIS camera zoom units range from 1 to 6000, which correspond to 1 to 6X zoom.
zoom.lm = lm(zoom.camera ~ zoom, data=data.frame(zoom=c(1,6000),
zoom.camera=c(1,6)))
summary(zoom.lm)
##
## Call:
## lm(formula = zoom.camera ~ zoom, data = data.frame(zoom = c(1,
## 6000), zoom.camera = c(1, 6)))
##
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9991665 NA NA NA
## zoom 0.0008335 NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 1 and 0 DF, p-value: NA
In this section we define models that are used to convert area and length between camera zoom levels to a common scale.
Download data for a reference object imaged at different zoom levels.
#if (!file.exists('zoom_calibration_data.txt')) {
# download.file('http://files.figshare.com/1845355/zoom_calibration_data.txt',
# 'zoom_calibration_data.txt')
#}
Read zoom calibrartion data
z.data = read.table(file="zoom_calibration_data.txt", sep="\t", header=TRUE)
Calculate px per cm
z.data$px_cm = z.data$length_px / z.data$length_cm
Convert LemnaTec zoom units to camera zoom units
z.data$zoom.camera = predict(object = zoom.lm, newdata=z.data)
vis.data$zoom = vis.data$sv_zoom
vis.data$sv.zoom.camera = predict(object = zoom.lm, newdata=vis.data)
vis.data$zoom = vis.data$tv_zoom
vis.data$tv.zoom.camera = predict(object = zoom.lm, newdata=vis.data)
Fit a variety of regression models to relative object area by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC)
Non-linear regression (exponential)
area.coef = coef(nls(log(rel_area) ~ log(a * exp(b * zoom.camera)),
z.data, start = c(a = 1, b = 0.01)))
area.coef = data.frame(a=area.coef[1], b=area.coef[2])
area.nls = nls(rel_area ~ a * exp(b * zoom.camera),
data = z.data, start=c(a=area.coef$a, b=area.coef$b))
summary(area.nls)
##
## Formula: rel_area ~ a * exp(b * zoom.camera)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.34409 0.01634 21.06 <2e-16 ***
## b 0.88510 0.01206 73.40 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2738 on 33 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 1.513e-08
Non-linear regression (polynomial)
area.pol = lm(rel_area ~ zoom.camera + I(zoom.camera^2), z.data)
summary(area.pol)
##
## Call:
## lm(formula = rel_area ~ zoom.camera + I(zoom.camera^2), data = z.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7458 -0.4521 0.1182 0.3087 1.6483
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.85322 0.54151 8.962 3.08e-10 ***
## zoom.camera -5.11648 0.45312 -11.292 1.07e-12 ***
## I(zoom.camera^2) 1.73120 0.08478 20.419 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.496 on 32 degrees of freedom
## Multiple R-squared: 0.9892, Adjusted R-squared: 0.9885
## F-statistic: 1465 on 2 and 32 DF, p-value: < 2.2e-16
AIC
AIC(area.nls, area.pol)
## df AIC
## area.nls 3 12.58188
## area.pol 4 55.10572
The exponential model minimizes AIC. Plot exponential model.
nls.plot = ggplot(z.data, aes(x=zoom.camera, y=rel_area)) +
geom_point(size=2.5) +
scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
scale_y_continuous(lim=c(0,17),
"Reference object relative pixel area") +
stat_smooth(data=z.data, aes(x=zoom.camera, y=rel_area),
method="nls", se=FALSE, formula=y ~ a * exp(b * x),
start=c(a=area.coef$a, b=area.coef$b)) +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
nls.plot = nls.plot + labs(title='Figure 12A')
print(nls.plot)
Fit a variety of regression models to px/cm by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC).
Non-linear regression (exponential)
len.coef = coef(nls(log(px_cm) ~ log(a * exp(b * zoom.camera)),
z.data[z.data$camera == 'VIS SV',], start = c(a = 1, b = 0.01)))
len.coef = data.frame(a=len.coef[1], b=len.coef[2])
len.nls = nls(px_cm ~ a * exp(b * zoom.camera),
data = z.data[z.data$camera == 'VIS SV',],
start=c(a=len.coef$a, b=len.coef$b))
summary(len.nls)
##
## Formula: px_cm ~ a * exp(b * zoom.camera)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 16.580319 0.271317 61.11 <2e-16 ***
## b 0.449384 0.004676 96.10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.035 on 15 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 1.958e-06
Length zoom correction using a 2 order polynomial (px/cm given a zoom setting). Length correction only works for side-view images right now because scale in top-view images are affected by the plant lifter position in addition to zoom.
len.poly = lm(px_cm ~ zoom.camera + I(zoom.camera^2),
data=z.data[z.data$camera == 'VIS SV',])
summary(len.poly)
##
## Call:
## lm(formula = px_cm ~ zoom.camera + I(zoom.camera^2), data = z.data[z.data$camera ==
## "VIS SV", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75431 -0.44016 -0.09013 0.40661 1.28788
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.1697 0.9308 33.49 9.13e-15 ***
## zoom.camera -8.9629 0.7901 -11.34 1.92e-08 ***
## I(zoom.camera^2) 6.5778 0.1503 43.75 2.24e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5886 on 14 degrees of freedom
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9995
## F-statistic: 1.718e+04 on 2 and 14 DF, p-value: < 2.2e-16
AIC
AIC(len.nls, len.poly)
## df AIC
## len.nls 3 53.28426
## len.poly 4 34.92045
The polynomial model minimizes AIC. Plot polynomial model.
poly.plot = ggplot(z.data[z.data$camera == 'VIS SV',], aes(x=zoom.camera, y=px_cm)) +
geom_point(size=2.5) +
scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
scale_y_continuous(lim=c(0,150),
"Reference object length scale (px/cm)") +
stat_smooth(data=z.data[z.data$camera == 'VIS SV',],
aes(x=zoom.camera, y=px_cm), method="lm",
formula=y ~ x + I(x^2)) +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
poly.plot = poly.plot + labs(title='Figure 12B')
print(poly.plot)
Create zoom-corrected VIS data frame
vis.data.zoom = vis.data[,c('plant_id', 'datetime', 'treatment', 'genotype', 'group', 'date', 'dap', 'solidity', 'outlier')]
Predict relative area zoom correction factors
vis.data$zoom.camera = vis.data$sv.zoom.camera
vis.data$sv_rel_area = predict(object = area.nls, newdata = vis.data)
vis.data$zoom.camera = vis.data$tv.zoom.camera
vis.data$tv_rel_area = predict(object = area.nls, newdata = vis.data)
Calculate total zoom-corrected side-view and top-view area
vis.data.zoom$sv_area = (vis.data$sv0_area / vis.data$sv_rel_area) + (vis.data$sv90_area / vis.data$sv_rel_area) +
(vis.data$sv180_area / vis.data$sv_rel_area) + (vis.data$sv270_area / vis.data$sv_rel_area)
vis.data.zoom$tv_area = vis.data$tv_area / vis.data$tv_rel_area
Calculate zoom-corrected lengths
vis.data$zoom.camera = vis.data$sv.zoom.camera
vis.data$px_cm = predict(object = len.poly, newdata=vis.data)
vis.data.zoom$extent_x = vis.data$extent_x / vis.data$px_cm
vis.data.zoom$extent_y = vis.data$extent_y / vis.data$px_cm
vis.data.zoom$height_above_bound = vis.data$height_above_bound / vis.data$px_cm
In this section we measure how well PlantCV estimates plant height.
Download data manually measured plant height data set (n = 173).
#if (!file.exists('height_model_data.txt')) {
# download.file('http://files.figshare.com/1845361/height_model_data.txt',
# 'height_model_data.txt')
#}
Read height data.
manual.ht.data = read.table(file="manual_height_samples.csv",sep=",",header=TRUE)
Convert human-readable date to Unix time
manual.ht.data$datetime = as.numeric(ymd_hms(manual.ht.data$timestamp, tz="America/Chicago"))
Get height data from the VIS data for each manual height sample
ht.data = merge(manual.ht.data, vis.data.zoom, by = c('plant_id', 'datetime'))
Convert LemnaTec zoom units to camera zoom units
ht.data$zoom.camera = predict(object = zoom.lm, newdata=ht.data)
Height linear model
height.model = lm(manual_height_cm~height_above_bound, ht.data)
summary(height.model)
##
## Call:
## lm(formula = manual_height_cm ~ height_above_bound, data = ht.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.99465 -0.23346 0.01155 0.24932 2.21212
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.372426 0.080998 -4.598 7.7e-06 ***
## height_above_bound 1.024824 0.002947 347.694 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6406 on 193 degrees of freedom
## Multiple R-squared: 0.9984, Adjusted R-squared: 0.9984
## F-statistic: 1.209e+05 on 1 and 193 DF, p-value: < 2.2e-16
Plot the height linear model.
hlm.plot = ggplot(ht.data, aes(x=height_above_bound, y=manual_height_cm)) +
geom_point(aes(color=factor(round(zoom.camera,1))),size=2.5) +
geom_smooth(method="lm", formula=y~x, color='black') +
scale_x_continuous(lim=c(0,60), "Estimated height (cm)") +
scale_y_continuous(lim=c(0,60), "Manually measured height (cm)") +
theme_bw() +
theme(legend.position=c(0.2,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Zoom")
hlm.plot = hlm.plot + labs(title='Figure 3A')
print(hlm.plot)
In this section we model fresh- and dry-weight above ground biomass using image measurements.
Download data manually measured plant biomass data set (n = 41).
#if (!file.exists('biomass_model_data.txt')) {
# download.file('http://files.figshare.com/1845360/biomass_model_data.txt',
# 'biomass_model_data.txt')
#}
Read biomass data.
manual.st.data = read.table(file='manual_biomass_samples.csv', sep=",", header=TRUE, stringsAsFactors=FALSE)
Get data from the VIS data for each manual biomass sample
st.data = merge(manual.st.data, vis.data.zoom, by = c('plant_id', 'datetime'))
Create out-of-frame indicator variable.
st.data$outind = NA
st.data[st.data$outlier == 'True',]$outind = 1
st.data[st.data$outlier == 'False',]$outind = 0
Genotype indicator variable.
st.data$group = NA
st.data[st.data$genotype == 'A10',]$group = 0
st.data[st.data$genotype == 'B100',]$group = 1
Full model. Includes side-view area, top-view area and height.
fw.full = lm(fresh_weight ~ sv_area * tv_area * height_above_bound, st.data)
Step-wise model selection with AIC.
fw.step = stepAIC(fw.full, direction="both")
## Start: AIC=-54.68
## fresh_weight ~ sv_area * tv_area * height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:tv_area:height_above_bound 1 0.031061 7.3442 -56.506
## <none> 7.3132 -54.680
##
## Step: AIC=-56.51
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area +
## sv_area:height_above_bound + tv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:tv_area 1 0.007503 7.3517 -58.464
## - tv_area:height_above_bound 1 0.059093 7.4033 -58.177
## - sv_area:height_above_bound 1 0.143762 7.4880 -57.711
## <none> 7.3442 -56.506
## + sv_area:tv_area:height_above_bound 1 0.031061 7.3132 -54.680
##
## Step: AIC=-58.46
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound +
## tv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - tv_area:height_above_bound 1 0.23599 7.5877 -59.169
## - sv_area:height_above_bound 1 0.31978 7.6715 -58.718
## <none> 7.3517 -58.464
## + sv_area:tv_area 1 0.00750 7.3442 -56.506
##
## Step: AIC=-59.17
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:height_above_bound 1 0.30538 7.8931 -59.551
## <none> 7.5877 -59.169
## + tv_area:height_above_bound 1 0.23599 7.3517 -58.464
## - tv_area 1 0.56816 8.1559 -58.208
## + sv_area:tv_area 1 0.18440 7.4033 -58.177
##
## Step: AIC=-59.55
## fresh_weight ~ sv_area + tv_area + height_above_bound
##
## Df Sum of Sq RSS AIC
## - tv_area 1 0.2631 8.156 -60.207
## <none> 7.893 -59.551
## + sv_area:height_above_bound 1 0.3054 7.588 -59.169
## + tv_area:height_above_bound 1 0.2216 7.671 -58.718
## + sv_area:tv_area 1 0.0478 7.845 -57.800
## - height_above_bound 1 0.8110 8.704 -57.541
## - sv_area 1 23.8156 31.709 -4.536
##
## Step: AIC=-60.21
## fresh_weight ~ sv_area + height_above_bound
##
## Df Sum of Sq RSS AIC
## <none> 8.156 -60.207
## + tv_area 1 0.263 7.893 -59.551
## - height_above_bound 1 0.654 8.810 -59.046
## + sv_area:height_above_bound 1 0.000 8.156 -58.208
## - sv_area 1 157.554 165.710 61.263
summary(fw.step)
##
## Call:
## lm(formula = fresh_weight ~ sv_area + height_above_bound, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.95619 -0.19476 -0.00613 0.06206 1.49454
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.455e-02 2.002e-01 0.173 0.864
## sv_area 3.972e-05 1.466e-06 27.093 <2e-16 ***
## height_above_bound -4.029e-02 2.309e-02 -1.745 0.089 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4633 on 38 degrees of freedom
## Multiple R-squared: 0.9842, Adjusted R-squared: 0.9833
## F-statistic: 1181 on 2 and 38 DF, p-value: < 2.2e-16
AIC model
fw.aic = lm(fresh_weight ~ sv_area + tv_area + height_above_bound +
sv_area*height_above_bound, st.data)
summary(fw.aic)
##
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound +
## sv_area * height_above_bound, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.14016 -0.10203 -0.00920 0.07919 1.56414
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.812e-01 2.356e-01 0.769 0.4468
## sv_area 4.538e-05 4.288e-06 10.584 1.34e-12 ***
## tv_area -1.026e-05 6.248e-06 -1.642 0.1093
## height_above_bound -6.232e-02 2.709e-02 -2.301 0.0273 *
## sv_area:height_above_bound 1.924e-07 1.598e-07 1.204 0.2366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4591 on 36 degrees of freedom
## Multiple R-squared: 0.9853, Adjusted R-squared: 0.9836
## F-statistic: 602 on 4 and 36 DF, p-value: < 2.2e-16
The AIC model contains tv_area and height which does not have a significant coefficient, test dropping.
fw.red = lm(fresh_weight ~ sv_area, st.data)
summary(fw.red)
##
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10964 -0.16717 0.00494 0.24429 1.30677
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.704e-01 1.002e-01 -2.70 0.0102 *
## sv_area 3.755e-05 7.931e-07 47.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4753 on 39 degrees of freedom
## Multiple R-squared: 0.9829, Adjusted R-squared: 0.9825
## F-statistic: 2241 on 1 and 39 DF, p-value: < 2.2e-16
Goodness of fit.
anova(fw.aic, fw.red)
## Analysis of Variance Table
##
## Model 1: fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area *
## height_above_bound
## Model 2: fresh_weight ~ sv_area
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 36 7.5877
## 2 39 8.8099 -3 -1.2222 1.9329 0.1417
SV area model.
sv.model = lm(fresh_weight ~ sv_area, st.data)
summary(sv.model)
##
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10964 -0.16717 0.00494 0.24429 1.30677
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.704e-01 1.002e-01 -2.70 0.0102 *
## sv_area 3.755e-05 7.931e-07 47.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4753 on 39 degrees of freedom
## Multiple R-squared: 0.9829, Adjusted R-squared: 0.9825
## F-statistic: 2241 on 1 and 39 DF, p-value: < 2.2e-16
Plot SV model
sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=fresh_weight)) +
geom_smooth(method="lm", color="black", formula = y ~ x) +
geom_point(size=2.5) +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
sv.model.plot = sv.model.plot + labs(title='Figure 4A')
print(sv.model.plot)
SV area model with out-of-frame.
sv.ind.model = lm(fresh_weight ~ sv_area + outind, st.data)
anova(sv.model, sv.ind.model)
## Analysis of Variance Table
##
## Model 1: fresh_weight ~ sv_area
## Model 2: fresh_weight ~ sv_area + outind
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 39 8.8099
## 2 38 7.0902 1 1.7196 9.2164 0.004316 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(sv.ind.model)
##
## Call:
## lm(formula = fresh_weight ~ sv_area + outind, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.78369 -0.19699 -0.01344 0.21652 1.20823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.202e-01 9.252e-02 -2.381 0.02241 *
## sv_area 3.831e-05 7.636e-07 50.174 < 2e-16 ***
## outind -5.242e-01 1.727e-01 -3.036 0.00432 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.432 on 38 degrees of freedom
## Multiple R-squared: 0.9862, Adjusted R-squared: 0.9855
## F-statistic: 1361 on 2 and 38 DF, p-value: < 2.2e-16
Plot SV model with out-of-frame.
sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight,
group=outlier, color=outlier)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Out-of-frame")
print(sv.model.out.plot)
SV area model with genotype
sv.gt.model = lm(fresh_weight ~ sv_area + group, st.data)
summary(sv.gt.model)
##
## Call:
## lm(formula = fresh_weight ~ sv_area + group, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24908 -0.29687 0.01704 0.21582 1.17723
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.150e-01 1.266e-01 -3.279 0.00223 **
## sv_area 3.775e-05 7.797e-07 48.414 < 2e-16 ***
## group 2.614e-01 1.460e-01 1.790 0.08135 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4624 on 38 degrees of freedom
## Multiple R-squared: 0.9842, Adjusted R-squared: 0.9834
## F-statistic: 1186 on 2 and 38 DF, p-value: < 2.2e-16
Plot SV model with genotype.
sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight,
group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype")
print(sv.model.gen.plot)
SV area model
dry.sv.model = lm(dry_weight ~ sv_area, st.data)
summary(dry.sv.model)
##
## Call:
## lm(formula = dry_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16009 -0.04743 0.01569 0.04033 0.18561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.915e-02 1.389e-02 -3.539 0.00106 **
## sv_area 4.425e-06 1.100e-07 40.233 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06591 on 39 degrees of freedom
## Multiple R-squared: 0.9765, Adjusted R-squared: 0.9759
## F-statistic: 1619 on 1 and 39 DF, p-value: < 2.2e-16
Plot dry-weight side-view model
dry.sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=dry_weight)) +
geom_smooth(method="lm", color="black", formula = y ~ x) +
geom_point(size=2.5) +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
dry.sv.model.plot = dry.sv.model.plot + labs(title='Supplemental Figure S2')
print(dry.sv.model.plot)
SV area model with out-of-frame.
dry.sv.ind.model = lm(dry_weight ~ sv_area + outind, st.data)
summary(dry.sv.ind.model)
##
## Call:
## lm(formula = dry_weight ~ sv_area + outind, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.112715 -0.030983 0.006637 0.024625 0.127446
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.052e-02 1.196e-02 -3.388 0.001650 **
## sv_area 4.556e-06 9.869e-08 46.166 < 2e-16 ***
## outind -9.023e-02 2.232e-02 -4.043 0.000248 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05583 on 38 degrees of freedom
## Multiple R-squared: 0.9836, Adjusted R-squared: 0.9827
## F-statistic: 1136 on 2 and 38 DF, p-value: < 2.2e-16
dry.sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight,
group=outlier, color=outlier)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Out-of-frame")
print(dry.sv.model.out.plot)
SV area model with genotypes
dry.sv.gt.model = lm(dry_weight ~ sv_area + group, st.data)
summary(dry.sv.gt.model)
##
## Call:
## lm(formula = dry_weight ~ sv_area + group, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15945 -0.04653 0.01473 0.03951 0.18609
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.015e-02 1.827e-02 -2.745 0.0092 **
## sv_area 4.426e-06 1.126e-07 39.315 <2e-16 ***
## group 1.813e-03 2.108e-02 0.086 0.9319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06676 on 38 degrees of freedom
## Multiple R-squared: 0.9765, Adjusted R-squared: 0.9752
## F-statistic: 788.8 on 2 and 38 DF, p-value: < 2.2e-16
dry.sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight,
group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype")
print(dry.sv.model.gen.plot)
In this section we analyze geometric traits from the complete VIS data set.
Remove plants that were sampled for biomass measurements.
vis.data.zoom = vis.data.zoom[!vis.data.zoom$plant_id %in% st.data$plant_id,]
Predict fresh- and dry-weight biomass from linear models
vis.data.zoom$fw_biomass = predict.lm(object = sv.model, newdata=vis.data.zoom)
vis.data.zoom$dw_biomass = predict.lm(object = dry.sv.model, newdata=vis.data.zoom)
Plot height for S. viridis and S. italica water treatments 100% and 33% full-capacity.
height.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
vis.data.zoom$genotype == 'B100') &
(vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33),],
aes(x=dap, y=height_above_bound, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,60),
name="Estimated height (cm)") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")
height.plot = height.plot + labs(title='Figure 3B')
print(height.plot)
Plot height for all genotypes.
height.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33,],
aes(x=dap, y=height_above_bound, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,60),
name="Estimated height (cm)") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)
height.facet.plot = height.facet.plot + labs(title='Supplemental Figure S1')
print(height.facet.plot)
For each day, calculate the mean and 95% confidence intervals for each genotype in the control water group.
height_per_day = function(vis.data) {
dap = c()
genotype = c()
int.low = c()
int.up = c()
est = c()
# Loop through each day and genotype
for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
if(d != 24) {
for(g in levels(factor(vis.data$genotype))) {
h.data = vis.data[vis.data$genotype == g & vis.data$treatment == 100 &
as.integer(vis.data$dap) == d,]$height_above_bound
if(sd(h.data) != 0) {
dap = c(dap,d)
genotype = c(genotype,g)
test = t.test(h.data)
int.low = c(int.low,test$conf.int[1])
int.up = c(int.up,test$conf.int[2])
est = c(est,test$estimate)
}
}
}
}
#drought_resp = drought_resp[-14,]
results = data.frame(dap=as.numeric(dap),
genotype=genotype,
conf.int.low=int.low,
conf.int.up=int.up,
mean=est)
return(results)
}
Height per day 95% CI. Write results to file.
height.perday.results = height_per_day(vis.data.zoom)
write.csv(height.perday.results,file="height.perday.100water.csv")
print(height.perday.results)
## dap genotype conf.int.low conf.int.up mean
## 1 11 A10 4.1138046 5.278839 4.696322
## 2 11 B100 2.5259319 3.326661 2.926297
## 3 11 R102 1.2684481 3.282102 2.275275
## 4 11 R128 1.6661243 2.966349 2.316237
## 5 11 R133 2.8262168 5.442689 4.134453
## 6 11 R161 2.8964352 3.279606 3.088021
## 7 11 R187 2.2374378 4.006971 3.122204
## 8 11 R20 4.2648901 6.073519 5.169204
## 9 11 R70 3.3847556 4.309511 3.847133
## 10 11 R98 5.5427726 6.565899 6.054336
## 11 12 A10 4.6895754 6.770786 5.730181
## 12 12 B100 2.7955379 4.609936 3.702737
## 13 12 R102 2.0129905 5.055558 3.534274
## 14 12 R128 2.3173878 4.551756 3.434572
## 15 12 R133 0.1833279 16.040642 8.111985
## 16 12 R161 2.7075759 5.056310 3.881943
## 17 12 R187 4.3345350 6.973768 5.654151
## 18 12 R20 -12.0782382 28.023729 7.972745
## 19 12 R70 3.0529634 5.786161 4.419562
## 20 12 R98 7.0338485 9.288841 8.161345
## 21 13 A10 7.2633774 8.933993 8.098685
## 22 13 B100 5.4971969 6.483222 5.990209
## 23 13 R102 3.9400526 6.147069 5.043561
## 24 13 R128 4.8671575 6.363675 5.615416
## 25 13 R133 7.5927485 9.939040 8.765894
## 26 13 R161 5.6093371 7.035664 6.322501
## 27 13 R187 5.3732471 7.338795 6.356021
## 28 13 R20 7.2176543 10.572329 8.894992
## 29 13 R70 6.2027656 7.171094 6.686930
## 30 13 R98 8.7080117 11.816909 10.262460
## 31 14 A10 7.2992964 10.650554 8.974925
## 32 14 B100 7.2095969 8.594076 7.901836
## 33 14 R102 4.6619405 8.705043 6.683492
## 34 14 R128 5.6149445 9.147871 7.381408
## 35 14 R133 11.8870482 12.800615 12.343831
## 36 14 R161 5.5132820 8.133326 6.823304
## 37 14 R187 7.7707801 10.776425 9.273602
## 38 14 R20 4.4397218 17.446651 10.943186
## 39 14 R70 6.2590147 10.235183 8.247099
## 40 14 R98 11.0626447 12.903724 11.983184
## 41 15 A10 10.5250697 12.024864 11.274967
## 42 15 B100 8.6034387 10.159156 9.381297
## 43 15 R102 5.5131064 8.935807 7.224457
## 44 15 R128 8.3515328 10.617965 9.484749
## 45 15 R133 13.4816282 16.926665 15.204147
## 46 15 R161 8.5636321 10.277182 9.420407
## 47 15 R187 7.8971054 11.444819 9.670962
## 48 15 R20 12.5540041 16.579996 14.567000
## 49 15 R70 9.4534249 11.399830 10.426628
## 50 15 R98 14.4840908 15.932854 15.208472
## 51 16 A10 10.3183950 13.664952 11.991674
## 52 16 B100 10.3217123 12.069260 11.195486
## 53 16 R102 7.4126348 10.542533 8.977584
## 54 16 R128 8.2558336 15.124237 11.690035
## 55 16 R133 16.4716388 17.833934 17.152786
## 56 16 R161 8.2040259 11.746160 9.975093
## 57 16 R187 11.1664495 13.991223 12.578836
## 58 16 R20 4.4901376 25.384677 14.937407
## 59 16 R70 10.1166241 15.935022 13.025823
## 60 16 R98 16.3400141 17.763702 17.051858
## 61 17 A10 13.8288699 16.273493 15.051181
## 62 17 B100 12.8021984 14.751262 13.776730
## 63 17 R102 4.6784026 20.775747 12.727075
## 64 17 R128 13.3720766 17.669471 15.520774
## 65 17 R133 17.4574601 17.473869 17.465665
## 66 17 R161 9.2299092 20.579302 14.904606
## 67 17 R187 13.3249948 15.732300 14.528647
## 68 17 R20 17.2492825 17.540747 17.395015
## 69 17 R70 13.8775205 16.001331 14.939426
## 70 17 R98 17.4389051 17.482331 17.460618
## 71 18 A10 14.4288846 16.848835 15.638860
## 72 18 B100 15.1035778 16.693924 15.898751
## 73 18 R102 11.9392515 17.839681 14.889466
## 74 18 R128 14.4065180 16.592639 15.499579
## 75 18 R161 11.3654767 16.590256 13.977866
## 76 18 R187 14.2918010 18.966135 16.628968
## 77 18 R20 17.3669600 17.559323 17.463141
## 78 18 R70 16.2359873 18.115003 17.175495
## 79 18 R98 17.4695785 17.486983 17.478281
## 80 19 A10 15.5274345 17.339111 16.433273
## 81 19 B100 14.1045628 16.056801 15.080682
## 82 19 R102 9.2740593 18.382349 13.828204
## 83 19 R128 17.4682044 17.481868 17.475036
## 84 19 R133 17.4691725 17.478305 17.473739
## 85 19 R161 13.7376679 17.711630 15.724649
## 86 19 R187 15.1000979 17.975410 16.537754
## 87 19 R20 17.4722809 17.482904 17.477592
## 88 19 R70 16.3674982 17.915366 17.141432
## 89 20 A10 15.4815153 17.904691 16.693103
## 90 20 B100 14.9940006 18.062583 16.528292
## 91 20 R102 14.8917746 18.071450 16.481612
## 92 20 R128 15.8643804 18.103082 16.983731
## 93 20 R133 17.4653321 17.482145 17.473739
## 94 20 R161 13.7669020 17.208033 15.487467
## 95 20 R187 16.8032817 17.759659 17.281470
## 96 20 R20 17.4540444 17.497470 17.475757
## 97 20 R70 15.5391179 18.384945 16.962032
## 98 20 R98 17.4675530 17.483962 17.475757
## 99 21 A10 18.0979070 22.866964 20.482435
## 100 21 B100 17.2659492 20.740679 19.003314
## 101 21 R102 10.3649359 27.239881 18.802408
## 102 21 R128 23.2259512 29.985433 26.605692
## 103 21 R133 35.2919365 38.324288 36.808112
## 104 21 R161 15.6427963 23.752576 19.697686
## 105 21 R187 20.9302871 28.030502 24.480395
## 106 21 R20 28.5799826 34.970515 31.775249
## 107 21 R70 22.9352632 26.500120 24.717692
## 108 21 R98 28.6623854 35.443911 32.053148
## 109 22 A10 20.8063040 25.919929 23.363117
## 110 22 B100 17.8425711 25.156676 21.499623
## 111 22 R102 16.7256650 28.193792 22.459728
## 112 22 R128 22.0557320 29.847211 25.951471
## 113 22 R133 33.3894751 49.007692 41.198583
## 114 22 R161 16.4579791 20.578905 18.518442
## 115 22 R187 25.3813857 35.138840 30.260113
## 116 22 R20 23.4710805 58.736805 41.103943
## 117 22 R70 19.0935859 36.220406 27.656996
## 118 22 R98 33.3321134 40.188754 36.760433
## 119 23 A10 23.7887258 29.592135 26.690430
## 120 23 B100 20.1451618 25.383980 22.764571
## 121 23 R102 11.4174069 32.037610 21.727508
## 122 23 R128 29.2037138 35.281459 32.242586
## 123 23 R133 41.3614576 45.130164 43.245811
## 124 23 R161 20.4242532 25.227423 22.825838
## 125 23 R187 23.1617592 38.261040 30.711400
## 126 23 R20 34.0321909 44.150975 39.091583
## 127 23 R70 30.2497912 34.894308 32.572050
## 128 23 R98 31.3077114 41.034772 36.171242
## 129 25 A10 28.8418538 34.316538 31.579196
## 130 25 B100 21.7374532 31.214009 26.475731
## 131 25 R102 19.3852591 38.594612 28.989935
## 132 25 R128 30.4325931 40.029886 35.231239
## 133 25 R133 47.2030660 54.937143 51.070105
## 134 25 R161 21.5395173 26.577798 24.058658
## 135 25 R187 35.8052978 41.428971 38.617134
## 136 25 R20 24.4029064 70.696036 47.549471
## 137 25 R70 28.2682884 48.061912 38.165100
## 138 25 R98 40.7567043 49.171868 44.964286
## 139 26 A10 31.3246339 38.418576 34.871605
## 140 26 B100 26.7993073 34.574345 30.686826
## 141 26 R102 14.9726330 46.784895 30.878764
## 142 26 R128 35.9731220 44.331988 40.152555
## 143 26 R133 52.6693551 52.719093 52.694224
## 144 26 R161 23.0900303 39.985494 31.537762
## 145 26 R187 33.1394226 43.142960 38.141191
## 146 26 R20 49.2489502 53.411285 51.330117
## 147 26 R70 41.7251538 46.738979 44.232066
## 148 26 R98 34.7437085 49.954721 42.349215
## 149 27 A10 35.1602826 41.245993 38.203138
## 150 27 B100 27.0240390 37.926865 32.475452
## 151 27 R102 23.9486786 47.091606 35.520142
## 152 27 R128 35.9592946 47.623371 41.791333
## 153 27 R133 52.6385536 52.932629 52.785592
## 154 27 R161 24.2442800 30.727013 27.485647
## 155 27 R187 41.2496516 44.943105 43.096378
## 156 27 R20 48.8256683 55.318931 52.072300
## 157 27 R70 34.9404743 52.806382 43.873428
## 158 27 R98 47.5987571 53.069042 50.333900
## 159 28 A10 37.3444960 43.747626 40.546061
## 160 28 B100 31.0852941 39.304128 35.194711
## 161 28 R102 19.7785402 51.925226 35.851883
## 162 28 R128 45.7837167 50.042465 47.913091
## 163 28 R133 52.6330116 52.823749 52.728380
## 164 28 R161 28.0109976 44.448898 36.229948
## 165 28 R187 35.2088327 47.003038 41.105935
## 166 28 R20 52.5767078 52.970566 52.773637
## 167 28 R70 49.9879362 52.206923 51.097429
## 168 28 R98 42.0724187 54.742865 48.407642
## 169 29 A10 40.2283979 45.252515 42.740457
## 170 29 B100 33.4612448 46.507641 39.984443
## 171 29 R102 26.5575297 53.548336 40.052933
## 172 29 R128 40.3678421 53.396169 46.882005
## 173 29 R133 52.6874775 53.051070 52.869274
## 174 29 R161 31.3730339 37.445685 34.409360
## 175 29 R187 44.0292757 47.520641 45.774958
## 176 29 R20 52.6909028 52.987872 52.839387
## 177 29 R70 42.0424398 56.607023 49.324731
## 178 29 R98 52.0544743 52.857213 52.455843
## 179 30 A10 41.1590692 48.054219 44.606644
## 180 30 B100 36.0961707 44.787331 40.441751
## 181 30 R102 27.1510111 54.451174 40.801093
## 182 30 R128 51.7339576 52.709507 52.221732
## 183 30 R133 52.6812412 52.877987 52.779614
## 184 30 R161 37.2649692 46.088567 41.676768
## 185 30 R187 39.0026350 49.389770 44.196202
## 186 30 R20 52.7731212 52.959449 52.866285
## 187 30 R70 52.6255605 52.942207 52.783884
## 188 30 R98 48.3607663 54.312277 51.336522
## 189 31 A10 44.9409479 48.995137 46.968042
## 190 31 B100 39.6784848 50.930006 45.304246
## 191 31 R102 32.3953484 54.156046 43.275697
## 192 31 R128 45.9896516 54.628337 50.308994
## 193 31 R133 52.8090833 53.001192 52.905138
## 194 31 R161 35.3843979 41.461346 38.422872
## 195 31 R187 46.3300941 49.725217 48.027656
## 196 31 R20 52.5361938 53.321900 52.929047
## 197 31 R70 45.1991354 56.510708 50.854922
## 198 31 R98 52.7570164 52.941683 52.849349
## 199 32 A10 43.9442034 51.699444 47.821823
## 200 32 B100 42.0255866 49.295184 45.660385
## 201 32 R102 33.9379036 55.911966 44.924935
## 202 32 R128 53.2778901 54.523848 53.900869
## 203 32 R133 54.3529911 54.469816 54.411403
## 204 32 R161 39.2900732 49.085036 44.187554
## 205 32 R187 43.5507291 49.570772 46.560751
## 206 32 R20 54.4236761 54.468606 54.446141
## 207 32 R70 54.3062174 54.417340 54.361779
## 208 32 R98 53.3227703 54.787243 54.055007
## 209 33 A10 48.7441481 52.620165 50.682157
## 210 33 B100 45.2839248 53.441351 49.362638
## 211 33 R102 38.4873097 55.935265 47.211287
## 212 33 R128 48.7761391 55.972918 52.374528
## 213 33 R133 54.3633798 54.459427 54.411403
## 214 33 R161 41.5367102 45.138371 43.337541
## 215 33 R187 47.6165089 50.414019 49.015264
## 216 33 R20 54.3071074 54.578858 54.442983
## 217 33 R70 46.8503717 57.968160 52.409266
## 218 33 R98 54.3418212 54.470459 54.406140
Statistical analysis for height differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.
analyze_height = function(vis.data, genotype) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(vis.data$dap)))) {
day = as.integer(day)
control = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 100,]$height_above_bound
drought = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 33,]$height_above_bound
test = t.test(x=control, y=drought)
days = c(days, day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}
Test for treatment effect on two-day intervals
a10.height.results = analyze_height(vis.data.zoom, 'A10')
b100.height.results = analyze_height(vis.data.zoom, 'B100')
Control for multiple testing by controlling the FDR
qvalues.height = p.adjust(c(a10.height.results$pvalue, b100.height.results$pvalue),method="fdr")
a10.height.results$qvalue = qvalues.height[1:nrow(a10.height.results)]
b100.height.results$qvalue = qvalues.height[
(nrow(a10.height.results)+1):(nrow(a10.height.results) + nrow(b100.height.results))]
Assign genotypes for merged table
a10.height.results$genotype = 'A10'
b100.height.results$genotype = 'B100'
print(a10.height.results)
## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 0.77273960 2.643169 9.582125e-04 0.0032431807 A10
## 2 12 -0.29914242 2.422404 1.206740e-01 0.2097086778 A10
## 3 13 -0.08842057 2.115835 7.039538e-02 0.1534676662 A10
## 4 14 -0.32917060 2.580711 1.239188e-01 0.2097086778 A10
## 5 15 0.28426142 2.919397 1.924394e-02 0.0529208463 A10
## 6 16 0.14666528 4.378048 3.748556e-02 0.0868086644 A10
## 7 17 -0.87910116 3.018279 2.505424e-01 0.3674622138 A10
## 8 18 0.48991495 3.897893 1.460564e-02 0.0428432116 A10
## 9 19 1.15428894 4.393747 1.873909e-03 0.0058894273 A10
## 10 20 2.29234606 6.796709 2.596630e-04 0.0009520977 A10
## 11 21 3.42568864 8.778005 6.793406e-05 0.0003736373 A10
## 12 22 5.05021107 10.781276 3.072841e-06 0.0001352050 A10
## 13 23 4.77919478 12.771985 2.155551e-04 0.0008622202 A10
## 14 25 3.75450813 10.565738 1.699507e-04 0.0008308701 A10
## 15 26 4.41239580 12.298773 1.911306e-04 0.0008409746 A10
## 16 27 5.94525469 13.501681 1.946187e-05 0.0001951542 A10
## 17 28 5.54257526 13.398589 5.606520e-05 0.0003736373 A10
## 18 29 5.90446577 12.852781 7.632409e-06 0.0001679130 A10
## 19 30 6.17854616 13.985667 2.217662e-05 0.0001951542 A10
## 20 31 6.60003839 14.577970 1.300773e-05 0.0001907800 A10
## 21 32 6.21299029 15.161443 6.320013e-05 0.0003736373 A10
## 22 33 1.73822539 19.572166 2.690657e-02 0.0696405324 A10
print(b100.height.results)
## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 -1.2594720 0.06992054 0.07673383 0.15346767 B100
## 2 12 -0.9886257 1.13640045 0.88658293 0.88658293 B100
## 3 13 -1.8408812 0.31467086 0.15621418 0.24547943 B100
## 4 14 -1.1123323 1.64762375 0.68794161 0.70394026 B100
## 5 15 -2.6455439 0.31836058 0.11604493 0.20970868 B100
## 6 16 -0.7864901 2.70009313 0.26085877 0.37025116 B100
## 7 17 -2.5824550 1.44347769 0.54836630 0.60320293 B100
## 8 18 -1.0746677 2.55761958 0.40054959 0.51835830 B100
## 9 19 -1.3873902 2.36907573 0.58999371 0.63316398 B100
## 10 20 -0.1953884 3.88187425 0.07458440 0.15346767 B100
## 11 21 -4.7055442 2.44925827 0.51436328 0.60320293 B100
## 12 22 -4.1699221 2.69196199 0.66036184 0.69180764 B100
## 13 23 -5.6805044 2.80840878 0.46303807 0.56593541 B100
## 14 25 -6.9139710 1.33553502 0.17774764 0.26968607 B100
## 15 26 -8.9304594 -0.52709508 0.02867192 0.07008691 B100
## 16 27 -7.3841917 1.12487918 0.14381616 0.23436707 B100
## 17 28 -8.7690652 0.63533080 0.08754828 0.16748367 B100
## 18 29 -6.1716464 2.70162092 0.43167803 0.54268095 B100
## 19 30 -6.3500065 2.33465144 0.35381611 0.48649716 B100
## 20 31 -2.6645326 4.95701544 0.54451691 0.60320293 B100
## 21 32 -2.8011143 5.16714436 0.54834969 0.60320293 B100
## 22 33 -3.1147374 7.53586129 0.38830762 0.51774350 B100
Output the A10 and B100 tables to the same file
write.table(a10.height.results, file='height.stats.csv', sep = ',',
row.names = FALSE, append = FALSE)
write.table(b100.height.results, file='height.stats.csv', sep = ',',
row.names = FALSE, append = TRUE, col.names = FALSE)
Plot fresh-weight biomass for S. viridis and S. italica water treatments 100% and 33% full-capacity.
biomass.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
vis.data.zoom$genotype == 'B100') &
(vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41),
name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")
biomass.plot = biomass.plot + labs(title='Figure 4B')
print(biomass.plot)
Plot fresh-weight biomass for all genotypes.
biomass.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33,],
aes(x=dap, y=fw_biomass, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41),
name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)
biomass.facet.plot = biomass.facet.plot + labs(title='Supplemental Figure S3')
print(biomass.facet.plot)
Statistical analysis for fresh-weight biomass differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.
analyze_biomass = function(vis.data, genotype) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(vis.data$dap)))) {
day = as.integer(day)
control = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 100,]$fw_biomass
drought = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 33,]$fw_biomass
# If the control group has a lot more replicates (e.g. A10 and B100), randomly sample to drought sample size
if (genotype == 'A10' | genotype == 'B100') {
control = sample(control, size = length(drought))
}
test = t.test(x=control, y=drought)
days = c(days, day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}
Test for treatment effect on two-day intervals
a10.biomass.results = analyze_biomass(vis.data.zoom, 'A10')
b100.biomass.results = analyze_biomass(vis.data.zoom, 'B100')
Control for multiple testing by controlling the FDR
qvalues.biomass = p.adjust(c(a10.biomass.results$pvalue, b100.biomass.results$pvalue),method="fdr")
a10.biomass.results$qvalue = qvalues.biomass[1:nrow(a10.biomass.results)]
b100.biomass.results$qvalue = qvalues.biomass[
(nrow(a10.biomass.results)+1):(nrow(a10.biomass.results) +
nrow(b100.biomass.results))]
Assign genotypes for merged table
a10.biomass.results$genotype = 'A10'
b100.biomass.results$genotype = 'B100'
print(a10.biomass.results)
## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 0.06618061 0.2366926 1.407098e-03 2.948206e-03 A10
## 2 12 -0.08602355 0.2214458 3.733164e-01 4.211774e-01 A10
## 3 13 -0.14878889 0.2530395 6.000685e-01 6.286432e-01 A10
## 4 14 -0.11355314 0.4706694 2.189056e-01 2.832896e-01 A10
## 5 15 -0.16532964 0.7004685 2.138495e-01 2.832896e-01 A10
## 6 16 -0.13466603 1.1361927 1.146885e-01 1.802249e-01 A10
## 7 17 -0.71306484 2.5798564 2.404504e-01 3.022805e-01 A10
## 8 18 1.01265567 2.7833441 1.984776e-04 5.822009e-04 A10
## 9 19 0.59986296 2.8235229 3.885584e-03 7.433292e-03 A10
## 10 20 1.51836098 4.7372845 6.640184e-04 1.623156e-03 A10
## 11 21 2.11154234 6.0531431 3.243150e-04 8.918661e-04 A10
## 12 22 4.28169719 8.8377380 9.412210e-06 3.185671e-05 A10
## 13 23 4.94273400 10.4089661 6.069847e-05 1.907666e-04 A10
## 14 25 6.47172203 11.7791467 6.402658e-07 3.130189e-06 A10
## 15 26 7.10409340 11.8459095 3.113334e-08 1.712334e-07 A10
## 16 27 8.85914954 12.3725610 6.213467e-12 4.556542e-11 A10
## 17 28 9.61069065 12.9472564 5.285230e-13 5.813753e-12 A10
## 18 29 8.91289326 12.3369483 3.464972e-12 3.049175e-11 A10
## 19 30 10.04685735 12.9896112 3.742783e-14 1.646825e-12 A10
## 20 31 10.14299516 13.5602211 3.490344e-13 5.119172e-12 A10
## 21 32 10.04460532 13.3093630 1.626586e-13 3.578490e-12 A10
## 22 33 9.97305515 14.9878564 1.469255e-06 6.464720e-06 A10
print(b100.biomass.results)
## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 -0.10273614 0.01603354 1.433897e-01 2.175569e-01 B100
## 2 12 -0.05331944 0.14763773 3.417760e-01 3.957407e-01 B100
## 3 13 -0.19433167 0.11252372 5.892521e-01 6.286432e-01 B100
## 4 14 -0.17223071 0.31125464 5.550045e-01 6.105049e-01 B100
## 5 15 -0.56141612 0.03872965 8.427560e-02 1.426203e-01 B100
## 6 16 -0.16839758 1.01001512 1.508052e-01 2.199402e-01 B100
## 7 17 -1.20761685 0.83516601 7.032150e-01 7.195688e-01 B100
## 8 18 -0.16852448 1.69244380 1.035239e-01 1.687056e-01 B100
## 9 19 -1.11131797 1.07162908 9.703722e-01 9.703722e-01 B100
## 10 20 0.95683425 3.18300525 8.908319e-04 2.062979e-03 B100
## 11 21 -1.02609382 3.09523737 3.076201e-01 3.759801e-01 B100
## 12 22 -0.71744932 4.15312089 1.549578e-01 2.199402e-01 B100
## 13 23 -1.72716481 7.54585618 1.698545e-01 2.335500e-01 B100
## 14 25 0.30686251 6.13591628 3.222730e-02 5.672006e-02 B100
## 15 26 -2.16265890 6.18863717 3.194544e-01 3.798917e-01 B100
## 16 27 4.43237220 9.24389988 8.819553e-06 3.185671e-05 B100
## 17 28 1.33865042 8.29025819 9.646199e-03 1.768470e-02 B100
## 18 29 2.22560476 8.60958217 2.396530e-03 4.793060e-03 B100
## 19 30 2.93818406 10.03007759 1.320144e-03 2.904316e-03 B100
## 20 31 8.43451895 12.61765594 5.634832e-10 3.541894e-09 B100
## 21 32 4.53402087 12.71198560 4.007013e-04 1.037109e-03 B100
## 22 33 8.36684530 14.24962413 1.865143e-06 7.460572e-06 B100
Output the A10 and B100 tables to the same file
write.table(a10.biomass.results, file='biomass.stats.csv', sep = ',',
row.names = FALSE, append = FALSE)
write.table(b100.biomass.results, file='biomass.stats.csv', sep = ',',
row.names = FALSE, append = TRUE, col.names = FALSE)
Calculate the difference in biomass per genotype per day between the 33% and 100% water treatment groups.
drought_response = function(vis.data, genotypes) {
# Initialize drought response data frame with days
drought_resp = data.frame(day=c(
min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))))
# Initialize genotypes
for(g in genotypes) {
control = paste(g,'control',sep='')
drought = paste(g,'drought',sep='')
# Calculate median biomass per treatment per day
drought_resp[,paste(control)] = 0
drought_resp[,paste(drought)] = 0
for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
drought_resp[drought_resp$day == d, paste(control)] =
median(vis.data[vis.data$genotype == g &
vis.data$treatment == 100 &
as.integer(vis.data$dap) == d,]$fw_biomass)
drought_resp[drought_resp$day == d, paste(drought)] =
median(vis.data[vis.data$genotype == g &
vis.data$treatment == 33 &
as.integer(vis.data$dap) == d,]$fw_biomass)
}
}
drought_resp = drought_resp[-14,]
# Calculate biomass loss to drought
days = c()
response = c()
genotype=c()
for(r in 1:nrow(drought_resp)) {
days = c(days, rep(drought_resp[r,]$day,length(genotypes)))
for(g in genotypes) {
control = paste(g,'control',sep='')
drought = paste(g,'drought',sep='')
response = c(response, drought_resp[r,paste(drought)] -
drought_resp[r,paste(control)])
genotype = c(genotype, g)
}
}
dr.df = data.frame(day=days,response=response,genotype=as.factor(genotype))
return(dr.df)
}
A10 vs B100 drought responses
drought.response.parents = drought_response(vis.data.zoom,c('A10','B100'))
Plot drought response results
genotype.colors = c('#E6A024','#4AB859')
resp.plot = ggplot(drought.response.parents,
aes(x=day, y=response, group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="loess") +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(name="Reduced accumulated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.8,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(values=genotype.colors) +
labs(color='Genotype')
resp.plot = resp.plot + labs(title='Figure 4C')
print(resp.plot)
In this experiment plants appear to follow relatively basic asymptotic growth patterns. Here we use non-linear least squares regression analysis to model a three-component logistic growth function for each genotype x treatment group.
First load three functions for non-linear growth curve analysis reproduced from:
Paine CET, Marthews TR, Vogt DR, Purves D, Rees M, Hector A, Turnbull LA (2012) How to fit nonlinear plant growth models and calculate growth rates: an update for ecologists. Methods in Ecology and Evolution 3: 245–256. doi: 10.1111/j.2041-210X.2011.00155.x.
output.logis.nlsList <- function(fit, times, CI = F, LOG = F, alpha = 0.05){
coef <- coef(fit)
params <- transform_param.logis(coef)
rates <- list()
groups <- rownames(params)
n.groups <- nrow(coef)
# compute rates for each group seperately
rates <- list()
for(i in 1:(n.groups)){
K <- params[i,1]; r <- params[i,2]; M0 <- params[i,3]
rates[[i]] = data.frame(
times = times,
M = (M0*K)/(M0+(K-M0)*exp(-r*times)),
AGR = (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
)
rates[[i]]$RGRt <- rates[[i]]$AGR/rates[[i]]$M
rates[[i]]$RGRm <- r*(1 - rates[[i]]$M/K)
if(LOG == T){
rates[[i]]$RGRt <- rates[[i]]$AGR
rates[[i]]$RGRm <- r*rates[[i]]$M*(1-rates[[i]]$M/K)
rates[[i]]$AGR <- rates[[i]]$AGR*exp(rates[[i]]$M)
}
# commute CIs for each group's estaimates, if desired
if(CI == T){
cov <- summary(fit)$cov[i,,]
x <- y <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=cov))
x$K <- y[,1]
x$r <- 1/y[,3]
x$M0 <- y[,1]/(1 + exp(y[,2]/y[,3]))
M <- AGR <- RGRt <- RGRm <- matrix(NA, ncol = length(times), nrow = nrow(x))
for(j in 1:nrow(x)){
K <- x[j,4]; r <- x[j,5]; M0 <- x[j,6]
M[j,] <- (M0*K)/(M0+(K-M0)*exp(-r*times))
AGR[j,] <- (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
RGRt[j,] <- AGR[j,]/M[j,]
RGRm[j,] <- r*(1 - M[j,]/K)
if(LOG ==T){
RGRt[j,] <- AGR[j,]
RGRm[j,] <- r*M[j,]*(1 - M[j,]/K)
AGR[j,] <- AGR[j,]*exp(M[j,])
}
}
}
CIs <- summarizer(list(M = M, AGR = AGR, RGRt = RGRt, RGRm = RGRm), alpha)
rates[[i]] <- cbind(rates[[i]], CIs)
}
names(rates) <- rownames(params)
# now compute differences among groups
diffs <- list()
for(i in 1:(n.groups-1)){
Ki <- params[i,1]; ri <- params[i,2]; M0i <- params[i,3]
for(j in (i+1):n.groups){
Kj <- params[j,1]; rj <- params[j,2]; M0j <- params[j,3]
diffs.ij = data.frame(
times = times,
diffM = rates[[i]]$M - rates[[j]]$M,
diffAGR = rates[[i]]$AGR - rates[[j]]$AGR,
diffRGRt = rates[[i]]$RGRt - rates[[j]]$RGRt
)
# comparing RGRm has to be done on a biomass basis. So it needs special treatment. First, we aev to know what range of biomasses are shared between groups.
Mmin <- max(min(rates[[i]]$M), min(rates[[j]]$M)) # yieds the range of overlapping masses between groups i & j
Mmax <- min(max(rates[[i]]$M), max(rates[[j]]$M))
#diffs.ij$Mseq <- seq(Mmin, Mmax, length = 100)
diffs.ij$Mseq <- seq(Mmin, Mmax, length = length(times))
if(LOG == F){
diffs.ij$diffRGRm <- ri*(1 - diffs.ij$Mseq/Ki) - rj*(1 - diffs.ij$Mseq/Kj)
} else{
diffs.ij$diffRGRm <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki) - rj*Mseq*(1 - diffs.ij$Mseq/Kj)
}
}
if(CI == T){
# get params for group i
covi <- summary(fit)$cov[i,,]
xi <- yi <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=covi))
xi$K <- yi[,1]
xi$r <- 1/yi[,3]
xi$M0 <- yi[,1]/(1 + exp(yi[,2]/yi[,3]))
# get params for group j
covj <- summary(fit)$cov[j,,]
xj <- yj <- data.frame(rmvnorm(n=1000, mean=c(coef[j, 1], coef[j, 2], coef[j, 3]), sigma=covj))
xj$K <- yj[,1]
xj$r <- 1/yj[,3]
xj$M0 <- yj[,1]/(1 + exp(yj[,2]/yj[,3]))
# now compute diffs for each random set of drawn parameters
Mi <- Mj <- AGRi <- AGRj <- RGRti <- RGRtj <- RGRmi <- RGRmj <- diffM <- diffAGR <- diffRGRt <- diffRGRm <- matrix(NA, ncol = length(times), nrow = nrow(xi))
for(k in 1:nrow(xi)){
Ki <- xi[k,4]; ri <- xi[k,5]; M0i <- xi[k,6]
Kj <- xj[k,4]; rj <- xj[k,5]; M0j <- xj[k,6]
Mi[k,] <- (M0i*Ki)/(M0i+(Ki-M0i)*exp(-ri*times))
Mj[k,] <- (M0j*Kj)/(M0j+(Kj-M0j)*exp(-rj*times))
AGRi[k,] <- (ri*M0i*Ki*(Ki-M0i)*exp(-ri*times))/(M0i+(Ki-M0i)*exp(-ri*times))^2
AGRj[k,] <- (rj*M0j*Kj*(Kj-M0j)*exp(-rj*times))/(M0j+(Kj-M0j)*exp(-rj*times))^2
RGRti[k,] <- AGRi[k,]/Mi[k,]
RGRtj[k,] <- AGRj[k,]/Mj[k,]
RGRmi[k,] <- ri*(1 - diffs.ij$Mseq/Ki)
RGRmj[k,] <- rj*(1 - diffs.ij$Mseq/Kj)
if(LOG == T){
RGRti[k,] <- AGRi[k,]
RGRtj[k,] <- AGRj[k,]
RGRmi[k,] <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki)
RGRmj[k,] <- rj*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Kj)
AGRi[k,] <- AGRi[k,]*exp(Mi[k,])
AGRj[k,] <- AGRj[k,]*exp(Mj[k,])
}
diffM[k,] <- Mi[k,] - Mj[k,]
diffAGR[k,] <- AGRi[k,] - AGRj[k,]
diffRGRt[k,] <- RGRti[k,] - RGRtj[k,]
diffRGRm[k,] <- RGRmi[k,] - RGRmj[k,]
}
CIs <- summarizer(list(diffM = diffM, diffAGR = diffAGR, diffRGRt = diffRGRt, diffRGRm = diffRGRm), alpha)
diffs[[paste(groups[i], groups[j], sep = "_")]] <- cbind(diffs.ij, CIs)
} else{
diffs[[paste(groups[i], groups[j], sep = "_")]] <- diffs.ij
}
} # end loop over pairwise combinations of groups
out <- list(params = params, rates = rates, diffs = diffs)
return(out)
}
transform_param.logis <- function(coef){
K = coef[1]
r = 1/(coef[3])
M0 = K/(1 + exp(coef[2]/coef[3])) #untransform best-fit parameters to K, r and M0
if(is.data.frame(K)){
out <- cbind(K, r, M0)
} else {
out <- c(K, r, M0)
}
names(out) <- c("K", "r", "M0")
return(out)
}
# this function returns confidence envelopes around growth trajectories, and growth rates.
summarizer <- function(dat, alpha){
n <- length(dat)
quantiles <- c(alpha/2, 1-(alpha/2))
CIs <- data.frame(matrix(NA, ncol(dat[[1]]), n*2))
names(CIs) <- paste(rep(names(dat), each = 2), c("lo", "hi"), sep = ".")
for(i in 1:n){
CIs[,(2*i-1):(2*i)] <- t(apply(dat[[i]], 2, quantile, quantiles, na.rm = T))
}
return(CIs)
}
Subset the VIS data
sub.vis.data = vis.data.zoom[(vis.data.zoom$treatment == 100 | vis.data.zoom$treatment == 33),
c("plant_id","dap","fw_biomass","group")]
sub.vis.data$group = factor(sub.vis.data$group)
Group data
grouped.sub.vis.data = groupedData(fw_biomass ~ dap | group, sub.vis.data)
Fit three-component logistic functions for each genotype x treatment group
fit.logis = nlsList(fw_biomass ~ SSlogis(dap, Asym, xmid, scal),
data = grouped.sub.vis.data)
Output estimated growth rate parameters at even time intervals
est.interval = seq(min(sub.vis.data$dap), max(sub.vis.data$dap), length = 100)
out.fit.logis = output.logis.nlsList(fit.logis, times = est.interval,
CI = TRUE, LOG = FALSE, alpha = 0.05)
Compute the time and magnitude of maximum growth rate
logis.results = data.frame(group=levels(sub.vis.data$group))
logis.results$max.AGR = 0
logis.results$max.AGR.dap = 0
logis.results$AGR.lo = 0
logis.results$AGR.hi = 0
group.order = names(out.fit.logis$rates)
for(i in 1:length(group.order)) {
logis.results[logis.results == group.order[i],]$max.AGR.dap =
out.fit.logis$rates[[i]]$times[out.fit.logis$rates[[i]]$AGR ==
max(out.fit.logis$rates[[i]]$AGR)]
logis.results[logis.results == group.order[i],]$max.AGR =
max(out.fit.logis$rates[[i]]$AGR)
logis.results[logis.results == group.order[i],]$AGR.lo =
out.fit.logis$rates[[i]]$AGR.lo[out.fit.logis$rates[[i]]$AGR ==
max(out.fit.logis$rates[[i]]$AGR)]
logis.results[logis.results == group.order[i],]$AGR.hi =
out.fit.logis$rates[[i]]$AGR.hi[out.fit.logis$rates[[i]]$AGR ==
max(out.fit.logis$rates[[i]]$AGR)]
}
print(logis.results)
## group max.AGR max.AGR.dap AGR.lo AGR.hi
## 1 A10-100 3.031821 24.23556 2.962031 3.099600
## 2 A10-33 1.800921 25.36054 1.724413 1.889975
## 3 B100-100 2.428901 26.48553 2.378344 2.479694
## 4 B100-33 1.693055 24.68555 1.604441 1.784634
## 5 R102-100 2.679992 28.06050 2.596513 2.757770
## 6 R102-33 1.890607 26.48553 1.813447 1.977527
## 7 R128-100 3.090714 25.58554 3.003102 3.174770
## 8 R128-33 1.847410 24.46055 1.770401 1.930073
## 9 R133-100 3.380220 23.11057 3.222844 3.526918
## 10 R133-33 2.292966 23.33557 2.180286 2.421985
## 11 R161-100 2.952080 25.58554 2.853000 3.043931
## 12 R161-33 1.813290 24.68555 1.729947 1.913824
## 13 R187-100 2.535501 24.46055 2.451728 2.629045
## 14 R187-33 1.676044 24.46055 1.596170 1.763532
## 15 R20-100 3.360632 22.88557 3.225079 3.499300
## 16 R20-33 1.949526 22.66058 1.858486 2.053219
## 17 R70-100 2.978971 24.01056 2.886229 3.083468
## 18 R70-33 1.654110 23.56057 1.571712 1.742591
## 19 R98-100 2.597831 24.23556 2.507745 2.691376
## 20 R98-33 1.538670 24.68555 1.467270 1.612429
Plot logistic growth models for S. viridis and S. italica and absolute growth rates over time
parent.groups = c("A10-100", "A10-33", "B100-100", "B100-33")
group.colors = c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")
# Biomass over time
gca.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
vis.data.zoom$genotype == 'B100') &
(vis.data.zoom$treatment == 100 | vis.data.zoom$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(group))) +
geom_point(size=2.5) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(name='Genotype-treatment',
values=group.colors, labels=parent.groups)
agr.plot = ggplot() +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(name="Absolute growth rate (g/day)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(name='Genotype-treatment',
values=group.colors, labels=parent.groups)
# Add logistic growth models and confidence intervals
for(g in 1:length(parent.groups)) {
group = parent.groups[g]
i = match(group, group.order)
group.rates = as.data.frame(out.fit.logis$rates[i])
col.name = gsub('-', '.', group)
conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
agr.conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.AGR.lo', sep='')],
rev(group.rates[,paste(col.name, '.AGR.hi', sep='')])))
gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.M', sep='')),
color=group.colors[g])
agr.plot = agr.plot +
geom_polygon(data=agr.conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.AGR', sep='')),
color=group.colors[g])
}
gca.plot = gca.plot + labs(title='Figure 4B')
print(gca.plot)
agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)
Plot logistic growth models for all Setaria genotypes
plot_gca = function(genotype, vis.data, out.fit.logis) {
groups = c(paste(genotype,'-100',sep=''), paste(genotype,'-33',sep=''))
group.colors = c("#00BFC4", "#F8766D")
# Biomass over time
gca.plot = ggplot(vis.data[(vis.data$genotype == genotype) &
(vis.data$treatment == 100 | vis.data$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(treatment))) +
geom_point(size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position='none',
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_grid(. ~ genotype)
# Add logistic growth models and confidence intervals
for(g in 1:length(groups)) {
group = groups[g]
i = match(group, group.order)
group.rates = as.data.frame(out.fit.logis$rates[i])
col.name = gsub('-', '.', group)
conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.M', sep='')),
color=group.colors[g])
}
return(gca.plot)
}
Multiple plot function modified from the Cookbook for R by Winston Chang: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
multiplot <- function(..., plotlist=NULL, cols=1, layout=NULL) {
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols),byrow=TRUE)
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
Generate the Setaria biomass/growth plots
a10.gca.plot = plot_gca('A10', vis.data.zoom, out.fit.logis)
b100.gca.plot = plot_gca('B100', vis.data.zoom, out.fit.logis)
r102.gca.plot = plot_gca('R102', vis.data.zoom, out.fit.logis)
r128.gca.plot = plot_gca('R128', vis.data.zoom, out.fit.logis)
r133.gca.plot = plot_gca('R133', vis.data.zoom, out.fit.logis)
r161.gca.plot = plot_gca('R161', vis.data.zoom, out.fit.logis)
r187.gca.plot = plot_gca('R187', vis.data.zoom, out.fit.logis)
r20.gca.plot = plot_gca('R20', vis.data.zoom, out.fit.logis)
r70.gca.plot = plot_gca('R70', vis.data.zoom, out.fit.logis)
r98.gca.plot = plot_gca('R98', vis.data.zoom, out.fit.logis)
Supplemental Figure S3
multiplot(a10.gca.plot, b100.gca.plot, r102.gca.plot, r128.gca.plot,
r133.gca.plot, r161.gca.plot, r187.gca.plot, r20.gca.plot,
r70.gca.plot, r98.gca.plot, cols=3)
agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)
In this section we combine data on the volume of water added to each plant per water application (up to the target weight) with fresh-weight biomass data to estimate water use efficiency.
Download water volume data (n = 33,496 water applications).
if (!file.exists('B2_watering_data_phenofront.csv')) {
download.file('http://files.figshare.com/1845363/B2_watering_data_phenofront.csv',
'B2_watering_data_phenofront.csv')
}
Read data and format for analysis
water.data = read.table(file="B2_watering_data_phenofront.csv", sep=",", header=TRUE)
Remove the empty pot controls (Barcodes start with zero).
water.data = water.data[grep('000A', water.data$plant.barcode, invert = TRUE),]
Remove snapshots that were not valid waterings (water.amount = -1).
water.data = water.data[water.data$water.amount != -1,]
Add a treatment group column coded in the barcode.
water.data$treatment <- NA
water.data$treatment[grep("AA", water.data$plant.barcode)] <- 100
water.data$treatment[grep("AB", water.data$plant.barcode)] <- 0
water.data$treatment[grep("AC", water.data$plant.barcode)] <- 16
water.data$treatment[grep("AD", water.data$plant.barcode)] <- 33
water.data$treatment[grep("AE", water.data$plant.barcode)] <- 66
Add a plant genotype column coded in the barcode.
water.data$genotype <- NA
water.data$genotype[grep("p1", water.data$plant.barcode)] <- 'A10'
water.data$genotype[grep("p2", water.data$plant.barcode)] <- 'B100'
water.data$genotype[grep("r1", water.data$plant.barcode)] <- 'R20'
water.data$genotype[grep("r2", water.data$plant.barcode)] <- 'R70'
water.data$genotype[grep("r3", water.data$plant.barcode)] <- 'R98'
water.data$genotype[grep("r4", water.data$plant.barcode)] <- 'R102'
water.data$genotype[grep("r5", water.data$plant.barcode)] <- 'R128'
water.data$genotype[grep("r6", water.data$plant.barcode)] <- 'R133'
water.data$genotype[grep("r7", water.data$plant.barcode)] <- 'R161'
water.data$genotype[grep("r8", water.data$plant.barcode)] <- 'R187'
Add a genotype x treatment group column.
water.data$group = paste(water.data$genotype,'-',water.data$treatment,sep='')
Encode the calendar time column as a date-time.
water.data$timestamp = ymd_hms(water.data$timestamp)
Add a column for days after planting since the planting date.
water.data$dap = as.numeric(water.data$timestamp - planting_date)
Some plants were removed from the system but were not inactivated, so remove them.
bad.cars = unique(water.data$car.tag[water.data$water.amount >120 &
water.data$dap > 14])
water.data = water.data[!water.data$car.tag %in% bad.cars,]
Extract water data for S. viridis for Figure 1B.
sv.water = water.data[(water.data$treatment == 100 | water.data$treatment == 33) &
water.data$genotype == 'A10',]
Classify water job as early or later in the day. Early water treatments started at ZT23, later treatments started at ZT8.
sv.water$early.late <- NA #restrict to the scheduled watering later in the run
sv.water$early.late[hour(sv.water$timestamp) < 7] <- 'ZT23'
sv.water$early.late[hour(sv.water$timestamp) < 15 &
hour(sv.water$timestamp) > 11 ] <- 'ZT8'
sv.water = sv.water[!is.na(sv.water$early.late),]
sv.water$dap = day(sv.water$timestamp) + 4
sv.water = sv.water[sv.water$dap > 14,]
Plot the S. viridis water volume data
water.plot = ggplot(sv.water, aes(y = water.amount, x = dap,
group = factor(interaction(early.late,
treatment,dap)),
fill = factor(interaction(early.late,
treatment)),
color = factor(interaction(early.late,
treatment)))) +
geom_boxplot(outlier.colour = NULL) +
scale_x_continuous("Days after planting") +
scale_y_continuous("Water volume (ml)") +
theme_bw() +
theme(legend.position = c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(fill='Water treatment', color='Water treatment')
water.plot = water.plot + labs(title='Figure 1B')
print(water.plot)
Calculate WUE for S. viridis and S. italica
wue.parents = function(water, biomass) {
# WUE data vectors
dap.list = c()
plant.list = c()
water.list = c()
biomass.list = c()
treatment.list = c()
group.list=c()
# Get unique barcodes for biomass
barcodes = unique(biomass[(biomass$genotype=='A10' | biomass$genotype=='B100') &
(biomass$treatment == 100 |
biomass$treatment == 33),]$plant_id)
for(barcode in barcodes) {
snapshots = biomass[biomass$plant_id==barcode,]
snapshots = snapshots[with(snapshots, order(dap)),]
for(row in 1:nrow(snapshots)) {
total.water = sum(water[water$dap <= snapshots[row,]$dap &
water$plant.barcode == barcode,]$water.amount)
dap.list = c(dap.list, snapshots[row,]$dap)
plant.list = c(plant.list, barcode)
water.list = c(water.list, total.water)
biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
treatment.list = c(treatment.list, snapshots[row,]$treatment)
group.list = c(group.list,snapshots[row,]$group)
}
}
wue.data = data.frame(plant_id=plant.list,
dap=dap.list,
water=water.list,
biomass=biomass.list,
treatment=treatment.list,
group=group.list)
return(wue.data)
}
parents.wue = wue.parents(water.data, vis.data.zoom)
parents.wue = parents.wue[parents.wue$water != 0,]
Plot S. viridis and S. italica WUE in mg/mL.
wue.plot = ggplot(parents.wue, aes(x=dap, y=(biomass/water)*1000,
color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess", size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-2.1, 26),
name="Water-use efficiency (mg/mL)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color='Genotype-treatment')
wue.plot = wue.plot + labs(title='Figure 4D')
print(wue.plot)
Genotype WUE function
wue = function(water, biomass, genotype) {
# WUE data vectors
dap.list = c()
plant.list = c()
water.list = c()
biomass.list = c()
treatment.list = c()
# Get unique barcodes for biomass
barcodes = unique(biomass[biomass$genotype==genotype &
(biomass$treatment == 100 |
biomass$treatment == 33),]$plant_id)
for(barcode in barcodes) {
snapshots = biomass[biomass$plant_id==barcode,]
snapshots = snapshots[with(snapshots, order(dap)),]
for(row in 1:nrow(snapshots)) {
total.water = sum(water[water$dap <= snapshots[row,]$dap &
water$plant.barcode == barcode,]$water.amount)
dap.list = c(dap.list, snapshots[row,]$dap)
plant.list = c(plant.list, barcode)
water.list = c(water.list, total.water)
biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
treatment.list = c(treatment.list, snapshots[row,]$treatment)
}
}
wue.data = data.frame(plant_id=plant.list,
dap=dap.list,
water=water.list,
biomass=biomass.list,
treatment=treatment.list)
return(wue.data)
}
Analyse pair-wise treatment differences for two-day increments for each genotype.
analyze_wue = function(wue.data) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(wue.data$dap)))) {
day = as.integer(day)
control = wue.data[(as.integer(wue.data$dap) == day |
as.integer(wue.data$dap) == day + 1) &
wue.data$treatment == 100,]
drought = wue.data[(as.integer(wue.data$dap) == day |
as.integer(wue.data$dap) == day + 1) &
wue.data$treatment == 33,]
control.wue = control$biomass / control$water
drought.wue = drought$biomass / drought$water
test = t.test(x=control.wue,y=drought.wue)
days = c(days,day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}
Analyze WUE for each genotype and aggregate results
wue.results = data.frame()
for(genotype in c('A10', 'B100')) {
genotype.wue = wue(water.data, vis.data.zoom, genotype)
genotype.wue = genotype.wue[genotype.wue$water != 0,]
genotype.wue.results = analyze_wue(genotype.wue)
genotype.wue.results$genotype = genotype
wue.results = rbind(wue.results, genotype.wue.results)
}
Control for multiple testing by controlling the FDR
qvalues.wue = p.adjust(wue.results$pvalue,method="fdr")
wue.results$qvalue = qvalues.wue
write.table(wue.results, file='wue.stats.csv', sep=',', row.names=FALSE)
print(wue.results)
## dap conf.int.low conf.int.up pvalue genotype qvalue
## 1 11 0.0001681857 7.003870e-04 0.0023589893 A10 0.03459851
## 2 12 -0.0002932779 6.855075e-04 0.4161121704 A10 0.79604067
## 3 13 -0.0003981046 7.310925e-04 0.5504252348 A10 0.85019455
## 4 14 -0.0003666556 1.086597e-03 0.3151858236 A10 0.65930472
## 5 15 -0.0006318583 1.290448e-03 0.4824013893 A10 0.85019455
## 6 16 -0.0011579704 1.505706e-03 0.7841660738 A10 0.85019455
## 7 17 -0.0030117768 2.033894e-03 0.6757645778 A10 0.85019455
## 8 18 -0.0019852804 1.537689e-03 0.7920861585 A10 0.85019455
## 9 19 -0.0022189870 1.423300e-03 0.6543284916 A10 0.85019455
## 10 20 -0.0018767476 2.487869e-03 0.7709120024 A10 0.85019455
## 11 21 -0.0011815956 3.333287e-03 0.3296523601 A10 0.65930472
## 12 22 -0.0011941099 3.447888e-03 0.3180544763 A10 0.65930472
## 13 23 -0.0030785137 4.736555e-03 0.6313454757 A10 0.85019455
## 14 25 -0.0017208810 2.529517e-03 0.6913194980 A10 0.85019455
## 15 26 -0.0016761719 2.370431e-03 0.7214693240 A10 0.85019455
## 16 27 -0.0015099451 1.887345e-03 0.8191576027 A10 0.85816511
## 17 28 -0.0016265463 1.712963e-03 0.9574064089 A10 0.95740641
## 18 29 -0.0016959640 1.418049e-03 0.8545599052 A10 0.87443339
## 19 30 -0.0017780029 1.240136e-03 0.7142953563 A10 0.85019455
## 20 31 -0.0018270439 1.050066e-03 0.5808623905 A10 0.85019455
## 21 32 -0.0018669052 9.419976e-04 0.5014284454 A10 0.85019455
## 22 33 -0.0019935110 2.547949e-03 0.7922267419 A10 0.85019455
## 23 11 -0.0004686284 1.653400e-04 0.3231119733 B100 0.65930472
## 24 12 -0.0003176706 4.659287e-04 0.6970305202 B100 0.85019455
## 25 13 -0.0006793007 2.052254e-04 0.2809421281 B100 0.65930472
## 26 14 -0.0005693316 7.964500e-04 0.7298612911 B100 0.85019455
## 27 15 -0.0013429980 2.693045e-04 0.1789099092 B100 0.56228829
## 28 16 -0.0011018089 1.726678e-03 0.6393283914 B100 0.85019455
## 29 17 -0.0053976734 -4.565038e-04 0.0249473645 B100 0.11089327
## 30 18 -0.0036398672 2.934236e-04 0.0894445265 B100 0.30273532
## 31 19 -0.0050384129 -6.422642e-04 0.0147704390 B100 0.08123741
## 32 20 -0.0044244661 -5.891914e-05 0.0447774175 B100 0.17910967
## 33 21 -0.0064021316 -1.321929e-03 0.0053896550 B100 0.03952414
## 34 22 -0.0057466915 -1.176396e-03 0.0051207769 B100 0.03952414
## 35 23 -0.0062955680 -1.663831e-03 0.0033048057 B100 0.03635286
## 36 25 -0.0058026119 -1.706030e-03 0.0009451266 B100 0.03381652
## 37 26 -0.0052016475 -1.378281e-03 0.0015371147 B100 0.03381652
## 38 27 -0.0043036325 -5.750897e-04 0.0125215176 B100 0.07870668
## 39 28 -0.0037893424 -2.728566e-04 0.0252030152 B100 0.11089327
## 40 29 -0.0033915507 1.520038e-04 0.0711932892 B100 0.26104206
## 41 30 -0.0027631015 6.920535e-04 0.2276453890 B100 0.65930472
## 42 31 -0.0025651583 7.463200e-04 0.2647434581 B100 0.65930472
## 43 32 -0.0025173692 8.428216e-04 0.3102535812 B100 0.65930472
## 44 33 -0.0029600383 2.256215e-03 0.7752593144 B100 0.85019455
In this section we use linear modeling to estimate tiller number.
Calculate height-width ratio.
vis.data.zoom$height_width_ratio = vis.data.zoom$height_above_bound / vis.data.zoom$extent_x
Download data for manually measured tiller counts 195 random images.
if (!file.exists('200tiller_counts.txt')) {
download.file('http://files.figshare.com/1845359/200tiller_counts.txt',
'200tiller_counts.txt')
}
Read manual tiller count data.
tiller.counts = read.table(file="200tiller_counts.txt",sep="\t",header=TRUE)
Format date.
tiller.counts$date = as.POSIXct(paste(paste(tiller.counts$year,
tiller.counts$month,
tiller.counts$day, sep='-'),
paste(tiller.counts$hours,
tiller.counts$minutes,
tiller.counts$seconds, sep=':'), sep=' '))
Merge with VIS snapshot table.
tiller.counts = merge(x=tiller.counts, y=vis.data.zoom, by = 'date')
Model tiller counts.
tiller.model.full = lm(ave_tillers ~ fw_biomass + height_above_bound + solidity +
extent_x + height_width_ratio, data=tiller.counts)
Variance inflation factor.
vif(tiller.model.full)
## fw_biomass height_above_bound solidity
## 16.615475 21.832350 1.591739
## extent_x height_width_ratio
## 14.738418 3.238483
Drop height and recalculate VIF
tiller.model1 = lm(ave_tillers ~ fw_biomass + solidity + extent_x +
height_width_ratio, data=tiller.counts)
vif(tiller.model1)
## fw_biomass solidity extent_x
## 10.014163 1.363769 11.118321
## height_width_ratio
## 1.847957
Drop extent x.
tiller.model2 = lm(ave_tillers ~ fw_biomass + solidity + height_width_ratio,
data=tiller.counts)
vif(tiller.model2)
## fw_biomass solidity height_width_ratio
## 1.024829 1.041759 1.039990
summary(tiller.model2)
##
## Call:
## lm(formula = ave_tillers ~ fw_biomass + solidity + height_width_ratio,
## data = tiller.counts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6877 -1.2636 -0.0169 0.8454 5.6053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.23276 0.54698 9.567 < 2e-16 ***
## fw_biomass 0.21998 0.01254 17.538 < 2e-16 ***
## solidity 0.15006 1.91657 0.078 0.938
## height_width_ratio -2.19370 0.32640 -6.721 2.12e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.707 on 187 degrees of freedom
## Multiple R-squared: 0.644, Adjusted R-squared: 0.6382
## F-statistic: 112.7 on 3 and 187 DF, p-value: < 2.2e-16
Drop solidity
tiller.model3 = lm(ave_tillers ~ fw_biomass + height_width_ratio, data=tiller.counts)
summary(tiller.model3)
##
## Call:
## lm(formula = ave_tillers ~ fw_biomass + height_width_ratio, data = tiller.counts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6886 -1.2680 -0.0274 0.8433 5.6022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.26031 0.41771 12.593 < 2e-16 ***
## fw_biomass 0.21986 0.01242 17.709 < 2e-16 ***
## height_width_ratio -2.18932 0.32072 -6.826 1.17e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.703 on 188 degrees of freedom
## Multiple R-squared: 0.6439, Adjusted R-squared: 0.6402
## F-statistic: 170 on 2 and 188 DF, p-value: < 2.2e-16
Plot partial regressions.
avPlots(model = tiller.model3, pch=16, layout=c(2,1), main="Figure 5B")
Download data for manually measured tiller counts for 646 random images.
if (!file.exists('660tiller_counts.txt')) {
download.file('http://files.figshare.com/1845358/660tiller_counts.txt',
'660tiller_counts.txt')
}
Predict tiller count for a second set of ~660 randomly selected plants.
tiller660 = read.table(file="660tiller_counts.txt", sep="\t", header=TRUE)
tiller660 = merge(x=tiller660, y=vis.data.zoom, by='datetime')
tiller660$predicted_tiller_count = predict.lm(object = tiller.model3,
newdata=tiller660)
Summarize the difference between the predicted and manual tiller counts
tiller.diff = tiller660$predicted_tiller_count - tiller660$tiller_count
t.test(tiller.diff)
##
## One Sample t-test
##
## data: tiller.diff
## t = 14.7453, df = 645, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.7978015 1.0429342
## sample estimates:
## mean of x
## 0.9203678
Plot manual versus predicted tiller counts.
tiller.plot = ggplot(data=tiller660, aes(x=tiller_count, y=predicted_tiller_count,
color=factor(genotype))) +
geom_point(size=2.5) +
geom_abline(intercept=0, slope=1) +
scale_x_continuous(lim=c(-1,14), "Manual tiller count") +
scale_y_continuous(lim=c(-1,14), "Predicted tiller count") +
theme_bw() +
theme(legend.position=c(0.45,0.95),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"),
legend.direction="horizontal") +
labs(color="Genotype")
tiller.plot = tiller.plot + labs(title='Figure 5C')
print(tiller.plot)
Predict tiller counts for the whole data set.
vis.data.zoom$tiller_count = predict.lm(object = tiller.model3, newdata = vis.data.zoom)
Plot tiller counts for S. viridis and S. italica water treatments 100% and 33% full-capacity.
tc.par.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
vis.data.zoom$genotype == 'B100') &
(vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33),],
aes(x=dap, y=tiller_count, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,12),
name="Estimated tiller count") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).
tc.par.plot = tc.par.plot + labs(title='Figure 5D')
print(tc.par.plot)
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).
Plot estimated tiller counts for all genotypes.
tiller.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33,],
aes(x=dap, y=tiller_count, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,12),
name="Estimated tiller count") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 3 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 4 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
tiller.facet.plot = tiller.facet.plot + labs(title='Supplemental Figure S4')
print(tiller.facet.plot)
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 3 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).
## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 4 rows containing missing values
## (geom_point).
## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
Output a table of all of the traits input and added here. Export all of the VIS traits (except color).
h.table = vis.data.zoom
h.table$plant_id = as.character(h.table$plant_id)
h.table$wue = NA
for(r in 1:nrow(h.table)) {
row = h.table[r,]
total.water = sum(water.data[water.data$dap <= row$dap & water.data$plant.barcode == row$plant_id,]$water.amount)
h.table[r,]$wue = row$sv_area / total.water
}
write.table(h.table,file="vis.traits.csv",quote=FALSE,sep=",",row.names=FALSE)